# All Hands on Deck: 
# The Role of Collaborative Platforms and Lead Organizations in Achieving Environmental Goals
# Heewon Lee and Yixin Liu 2024
# Journal of Public Administration Research and Theory
# Replication File

####Packages####
library(dplyr)
library(tidyverse)
library(data.table)
library(ggplot2)
library(cowplot)
library(egg)
library(gridExtra)
library(grid)
library(ggpubr)
library(maps)
library(usmap)
library(ggmap)
library(mapproj)
library(USAboundariesData)
library(USAboundaries)
library(housingData)
library(fixest)
library(ggiplot)
library(did2s)
library(xtable)

####Functions####
# ggplot titlebox
element_textbox <- function(...) {
  el <- element_text(...)
  class(el) <- c("element_textbox", class(el))
  el
}

element_grob.element_textbox <- function(element, ...) {
  text_grob <- NextMethod()
  rect_grob <- element_grob(calc_element("strip.background", theme_bw()))
  
  ggplot2:::absoluteGrob(
    grid::gList(
      element_grob(calc_element("strip.background", theme_bw())),
      text_grob
    ),
    height = grid::grobHeight(text_grob), 
    width = grid::unit(1, "npc")
  )
}

# Two decimal places
scaleFUN <- function(x) {
  return(sprintf("%.2f", x))
}

# generate plots
dv_colors = list(AQI = '#56B4E9', BDR = '#9999CC')
dvs = c("AQI", "BDR")

did2s_function <- function(data, dvs, dv_colors, ccc_type = "All", coeff_range = NULL, custom_title_part = NULL) {
  require(ggplot2)
  require(did)
  require(fixest)
  
  plots <- list() # Store plots
  overall_effects <- list() # Store overall effect regressions
  
  for(dv in dvs) {
    # Run regression for event plot
    es <- did2s(
      data = data,
      yname = dv, 
      first_stage = ~ pop_density + 
        white + 
        freight +
        laborforce +
        annual_temp +
        annual_prec +
        total_incen_no +
        hydro_generator_per +
        hydro_capacity_per +
        coal_generator_per +
        coal_capacity_per +
        Total_Generators +
        Total_Capacity +
        total_regulation +
        violation +
        pro_env
      | fips + year,
      second_stage = ~ i(rel_time, ref = c(0, Inf)), 
      treatment = "treated",
      cluster_var = "state"
    )
    
    # Obtain standard errors for event plot
    es_summary <- summary(es, conley(100))
    coeftable <- es_summary$coeftable
    coefficients <- coeftable[, "Estimate"]
    std_errors <- coeftable[, "Std. Error"]
    
    # Create a dataframe for event plot
    df_plot <- data.frame(
      term = as.numeric(sub("rel_time::", "", rownames(coeftable))),
      estimate = coefficients,
      conf_low = coefficients - 1.96 * std_errors,
      conf_high = coefficients + 1.96 * std_errors
    )
    
    # Filter coefficients for event plot if coeff_range is provided
    if (!is.null(coeff_range)) {
      df_plot <- df_plot[df_plot$term >= coeff_range[1] & df_plot$term <= coeff_range[2], ]
    }
    
    df_zero <- data.frame(term = 0, estimate = 0, conf_low = NA, conf_high = NA)
    df_plot <- rbind(df_zero, df_plot)
    
    # Construct the event plot title
    title_prefix <- ifelse(is.null(custom_title_part), '', paste0(custom_title_part, ' '))
    plot_title <- paste0(title_prefix, '\nEstimate Change of ', dv, ' and 95% CIs')
    
    # Plot event plot
    plot <- ggplot(df_plot, aes(x = term, y = estimate, ymin = conf_low, ymax = conf_high)) +
      geom_vline(aes(xintercept = 0), color = "#CC6666", linetype = "dashed", size = 1) +
      geom_hline(aes(yintercept = 0), color = "#CC6666", size = 1) +
      geom_point(aes(color = dv), size = 2) +
      geom_line(aes(color = dv, group = 1), linewidth = 0.7) + 
      geom_errorbar(aes(color = dv), width = 0.7) +
      scale_color_manual(values = dv_colors[[dv]]) +
      labs(
        title = plot_title,
        x = 'Years Relative to CCC Adoption',
        y = ''
      ) +
      theme_bw() +
      theme(
        plot.title = element_textbox(
          size = 16, hjust = 0.5, margin = margin(t = 5, b = 5)
        ),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.title.x = element_text(size = 12),
        legend.position = "none"
      ) +
      scale_x_continuous(limits = c(-15.5, 20.5), breaks = c(-15, -10, -5, 0, 5, 10, 15, 20))
    
    plots[[dv]] <- plot
    
    # Run overall effect regression
    overall_effect <- did2s(
      data = data,
      yname = dv, 
      first_stage = ~ pop_density + 
        white + 
        freight +
        laborforce +
        annual_temp +
        annual_prec +
        total_incen_no +
        hydro_generator_per +
        hydro_capacity_per +
        coal_generator_per +
        coal_capacity_per +
        Total_Generators +
        Total_Capacity +
        total_regulation +
        violation +
        pro_env
      | fips + year,
      second_stage = ~ i(treated, ref = FALSE),
      treatment = "treated",
      cluster_var = "state"
    )
    
    # Store the summary of overall effect regression
    overall_effects[[dv]] <- summary(overall_effect, conley(100))
  }
  
  return(list("Plots" = plots, "ATT" = overall_effects))
}

did2s_function_placebo <- function(data, dvs, dv_colors, ccc_type = "All", coeff_range = NULL, custom_title_part = NULL) {
  require(ggplot2)
  require(did)
  require(fixest)
  
  plots <- list() # Store plots
  
  for(dv in dvs) {
    # Run regression
    es <- did2s(
      data = data,
      yname = dv, 
      first_stage = ~ pop_density + 
        white + 
        freight +
        laborforce +
        annual_temp +
        annual_prec +
        total_incen_no +
        hydro_generator_per +
        hydro_capacity_per +
        coal_generator_per +
        coal_capacity_per +
        Total_Generators +
        Total_Capacity +
        total_regulation +
        violation +
        pro_env
      | fips + year,
      second_stage = ~ i(rel_time_placebo, ref = c(0, Inf)), 
      treatment = "treated_placebo",
      cluster_var = "state"
    )
    
    # Obtain standard errors
    es_summary <- summary(es, conley(100))
    coeftable <- es_summary$coeftable
    coefficients <- coeftable[, "Estimate"]
    std_errors <- coeftable[, "Std. Error"]
    
    # Create a dataframe for ggplot
    df_plot <- data.frame(
      term = as.numeric(sub("rel_time_placebo::", "", rownames(coeftable))),
      estimate = coefficients,
      conf_low = coefficients - 1.96 * std_errors,
      conf_high = coefficients + 1.96 * std_errors
    )
    
    # Filter coefficients if coeff_range is provided
    if (!is.null(coeff_range)) {
      df_plot <- df_plot[df_plot$term >= coeff_range[1] & df_plot$term <= coeff_range[2], ]
    }
    
    df_zero <- data.frame(term = 0, estimate = 0, conf_low = NA, conf_high = NA)
    df_plot <- rbind(df_zero, df_plot)
    
    # Construct the plot title
    title_prefix <- ifelse(is.null(custom_title_part), '', paste0(custom_title_part, ' '))
    plot_title <- paste0(title_prefix, '\nEstimate Change of ', dv, ' and 95% CIs')
    
    # Set default or custom x-axis limits and breaks
    if (is.null(coeff_range)) {
      x_limits <- c(-10.5, 2.5)
      x_breaks <- c(-10, -5, 0, 2)
    } else {
      x_limits <- c(min(coeff_range[1], -10.5), max(coeff_range[2], 2.5))
      x_breaks <- seq(from = min(coeff_range[1], -10), to = max(coeff_range[2], 2), by = 5)
    }
    
    # Plotting
    plot <- ggplot(df_plot, aes(x = term, y = estimate, ymin = conf_low, ymax = conf_high)) +
      geom_vline(aes(xintercept = 0), color = "#CC6666", linetype = "dashed", size = 1) +
      geom_hline(aes(yintercept = 0), color = "#CC6666", size = 1) +
      geom_point(aes(color = dv), size = 2) +
      geom_line(aes(color = dv, group = 1), linewidth = 0.7) + 
      geom_errorbar(aes(color = dv), width = 0.7) +
      scale_color_manual(values = dv_colors[[dv]]) +
      labs(
        title = plot_title,
        x = 'Years Relative to Placebo CCC Adoption',
        y = ''
      ) +
      theme_bw() +
      theme(
        plot.title = element_textbox(
          size = 16, hjust = 0.5, margin = margin(t = 5, b = 5)
        ),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.title.x = element_text(size = 12),
        legend.position = "none"
      ) +
      scale_x_continuous(limits = x_limits, breaks = x_breaks)
    
    plots[[dv]] <- plot
  }
  
  return(plots)
}

# descriptive table
data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,
                  varname)
  data_sum <- rename(data_sum, c("mean" = varname))
  return(data_sum)
}

####Data####
# Setup working directory
WD = getwd()
setwd(WD)

df<-read.csv("df_jpart.csv")
df_fullmap<-read.csv("df_fullmap_jpart.csv")

####Subgroup####
# lead organization subgroups
df_gov<-subset(df,df$orgtype!=4)
df_state<-subset(df,df$orgtype==1 | df$orgtype==10)
df_local<-subset(df,df$orgtype==2 | df$orgtype==10)
df_regional<-subset(df,df$orgtype==3 | df$orgtype==10)
df_npo<-subset(df,df$orgtype==4 | df$orgtype==10)

# mechanism subgroups
df_goal<-subset(df,df$goal==1 | df$orgtype==10)

df_resource_donation<-subset(df,df$resources_donation | df$orgtype==10)
df_resource_fedfund<-subset(df,df$resources_fedfund | df$orgtype==10)

df_group_crosssector<-subset(df,df$group_crosssector==1 | df$orgtype==10)
df_group_open<-subset(df,df$group_open==1 | df$orgtype==10)
df_group_benefit<-subset(df,df$group_member_benefit | df$orgtype==10)
df_group_resp<-subset(df,df$group_member_resp | df$orgtype==10)

####Map####
states <- plot_usmap(exclude = c("AK","HI","DC"), 
                     color = "yellow",
                     fill = alpha(0.01))

counties<- plot_usmap(data = df_fullmap,
                      exclude = c("AK","HI","DC"), 
                      values = "orgtype",
                      color = "black",
                      size = 0.1)

# map for all counties
counties[[1]]$orgtype <- factor(counties[[1]]$orgtype, levels = c("Never Treated",
                                                                  "Local-led",
                                                                  "Nonprofit-led",
                                                                  "Regional-led",
                                                                  "State-led"))

fullmap_plot<-ggplot() +  
  geom_polygon(data=counties[[1]], 
               aes(x=x, 
                   y=y, 
                   group=group, 
                   fill = counties[[1]]$orgtype), 
               color = "black",
               size = 0.1) +  
  geom_polygon(data=states[[1]], 
               aes(x=x, 
                   y=y, 
                   group=group), 
               color = "white", 
               fill = alpha(0.01)) + 
  labs(x=element_blank(), y=element_blank(), 
       title="U.S. Counties by CCC Treatment Conditions") +
  coord_equal() +
  theme_map() +
  theme_bw() +
  theme(legend.position=c(0.12, 0.15),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.text = element_text(size = 14),
        axis.text=element_blank(),
        axis.ticks = element_blank(),
        strip.background = element_blank(),
        strip.text.x = element_blank(),
        plot.title = element_textbox(
          size = 22, hjust = 0.5, margin = margin(t = 5, b = 5)))
fullmap_plot

# counties observed
df_observed_map<-df[,c(3,4,10,11,6)]
df_observed_map$fips<-sprintf("%05d", df_observed_map$fips)
df_observed_map$fips<-as.character(df_observed_map$fips)
df_observed_map$orgtype1[df_observed_map$orgtype==1]<-"State-led"
df_observed_map$orgtype1[df_observed_map$orgtype==2]<-"Local-led"
df_observed_map$orgtype1[df_observed_map$orgtype==3]<-"Regional-led"
df_observed_map$orgtype1[df_observed_map$orgtype==4]<-"Nonprofit-led"
df_observed_map$orgtype1[df_observed_map$orgtype==10]<-"Never Treated"

counties_observed<- plot_usmap(data = df_observed_map, 
                               exclude = c("AK","HI","DC"),
                               values = "orgtype1", 
                               color = "black",
                               size = 0.1)

counties_observed[[1]]$orgtype1 <- factor(counties_observed[[1]]$orgtype1, 
                                     levels = c("Never Treated",
                                                "Local-led",
                                                "Nonprofit-led",
                                                "Regional-led",
                                                "State-led"))

observedmap_plot<-ggplot() +  
  geom_polygon(data=counties_observed[[1]], 
               aes(x=x, 
                   y=y, 
                   group=group, 
                   fill = counties_observed[[1]]$orgtype1), 
               color = "black",
               size = 0.1) +  
  geom_polygon(data=states[[1]], 
               aes(x=x, 
                   y=y, 
                   group=group), 
               color = "white", 
               fill = alpha(0.01)) + 
  labs(x=element_blank(), y=element_blank(), 
       title="Sample Counties by CCC Treatment Conditions") +
  coord_equal() +
  theme_map() +
  theme_bw() +
  scale_fill_discrete(label=c("Never Treated","Local-led","Nonprofit-led",
                              "Regional-led","State-led",
                              "Air Quality Not Available")) +
  theme(legend.position=c(0.18, 0.18),
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.text = element_text(size = 14),
        axis.text=element_blank(),
        axis.ticks = element_blank(),
        strip.background = element_blank(),
        strip.text.x = element_blank(),
        plot.title = element_textbox(
          size = 22, hjust = 0.5, margin = margin(t = 5, b = 5)))
observedmap_plot

#####Figure 3####
figure3 <- ggarrange(fullmap_plot, observedmap_plot,ncol = 2, nrow = 1)
figure3<-plot_grid(figure3)
figure3
ggsave(file = "figure3.jpg", width = 15, height = 5.5, dpi = 666)
ggsave(file = "figure3.eps", width = 15, height = 5.5, dpi = 666)

####Adoption Pattern####
df_ccc<-df%>%
  dplyr::group_by(coalitionid) %>%
  dplyr::summarise(ccc=mean(adoption),
            orgtype=dplyr::first(orgtype))

df_ccc$ccc_text<-as.character(df_ccc$ccc)
df_ccc<-subset(df_ccc, df_ccc$ccc!=0)

df_ccc$orgtype[df_ccc$orgtype==1]<-"State-led"
df_ccc$orgtype[df_ccc$orgtype==2]<-"Local-led"
df_ccc$orgtype[df_ccc$orgtype==3]<-"Regional-led"
df_ccc$orgtype[df_ccc$orgtype==4]<-"Nonprofit-led"

#####Figure 4####
figure4 <- ggplot(df_ccc,
                   aes(x=ccc_text,
                       fill=orgtype)) + 
  geom_bar(color="black") + 
  labs(x=element_blank(), y=element_blank(),
       title="Number of CCC Adoptions (1990-2020)") +
  theme_bw() +
  theme(legend.title=element_blank(),
        legend.text = element_text(size=14),
        legend.background = element_blank(),
        axis.text=element_text(size=10),
        legend.justification=c(1,1), legend.position=c(1,1),
        plot.title = element_textbox(
          size = 14, hjust = 0.5, margin = margin(t = 5, b = 5))) +
  scale_fill_manual(values=c("#A3A500","#00BF7D","#00B0F6", "#E76BF3"))

figure4
ggsave(file = "figure4.jpg", width = 8, height = 4, dpi = 666)
ggsave(file = "figure4.eps", width = 8, height = 4, dpi = 666)

####Overall Effects####
plots_out <- did2s_function(df, dvs, dv_colors, custom_title_part = "All CCCs:")

#####Table 2####
mean(df$AQI[df$treated==0 & df$ccc==1])
mean(df$BDR[df$treated==0 & df$ccc==1])

aqi_att<-plots_out$ATT$AQI
aqi_att

bad_att<-plots_out$ATT$BDR
bad_att

etable(aqi_att, bad_att,
       fitstat=~n+ar2, digits =3, digits.stats = 3, tex = TRUE)

#####Figure 5####
aqi_plot<-plots_out$Plots$AQI + ggplot2::scale_y_continuous(limits=c(-6,2.5),labels=scaleFUN)
aqi_plot

bad_plot<-plots_out$Plots$BDR + ggplot2::scale_y_continuous(limits=c(-0.05,0.02),labels=scaleFUN)
bad_plot


figure5<-ggarrange(aqi_plot, bad_plot, 
                   ncol = 2, nrow = 1)

figure5<-plot_grid(figure5)
figure5
ggsave(file = "figure5.jpg", width = 12, height = 4, dpi = 666)
ggsave(file = "figure5.eps", width = 12, height = 4, dpi = 666)

####Nonprofit-led vs Government-led####
mean(df_gov$AQI[df_gov$treated==0 & df_gov$ccc==1])
mean(df_gov$BDR[df_gov$treated==0 & df_gov$ccc==1])
mean(df_npo$AQI[df_npo$treated==0 & df_npo$ccc==1])
mean(df_npo$BDR[df_npo$treated==0 & df_npo$ccc==1])

# Nonprofit-led CCCs
plots_out_npo <- did2s_function(df_npo, dvs, dv_colors, 
                              ccc_type = "Nonprofit-led",
                              custom_title_part = "Nonprofit-led CCCs:")

aqi_att_npo<-plots_out_npo$ATT$AQI
aqi_att_npo

bad_att_npo<-plots_out_npo$ATT$BDR
bad_att_npo

aqi_plot_npo<-plots_out_npo$Plots$AQI + ggplot2::scale_y_continuous(limits=c(-8,4),labels=scaleFUN)
aqi_plot_npo

bad_plot_npo<-plots_out_npo$Plots$BDR + ggplot2::scale_y_continuous(limits=c(-0.06,0.03),labels=scaleFUN)
bad_plot_npo

# Government-led CCCs
plots_out_gov <- did2s_function(df_gov, dvs, dv_colors, 
                              ccc_type = "Government-led",
                              custom_title_part = "Government-led CCCs:")

aqi_att_gov<-plots_out_gov$ATT$AQI
aqi_att_gov

bad_att_gov<-plots_out_gov$ATT$BDR
bad_att_gov

aqi_plot_gov<-plots_out_gov$Plots$AQI + ggplot2::scale_y_continuous(limits=c(-8,4),labels=scaleFUN)
aqi_plot_gov

bad_plot_gov<-plots_out_gov$Plots$BDR + ggplot2::scale_y_continuous(limits=c(-0.06,0.03),labels=scaleFUN)
bad_plot_gov

#####Table 3####
etable(aqi_att_gov,bad_att_gov, 
       aqi_att_npo,bad_att_npo,
       fitstat=~n+ar2, digits =3, digits.stats = 3, tex = TRUE)

#####Figure 6####
figure6<-ggarrange(aqi_plot_gov, bad_plot_gov,
                    aqi_plot_npo, bad_plot_npo,
                    ncol = 2, nrow = 2)

figure6<-plot_grid(figure6)
figure6
ggsave(file = "figure6.jpg", width = 12, height = 8, dpi = 666)
ggsave(file = "figure6.eps", width = 12, height = 8, dpi = 666)

####Government Subgroups####
mean(df_state$AQI[df_state$treated==0 & df_state$ccc==1])
mean(df_state$BDR[df_state$treated==0 & df_state$ccc==1])
mean(df_regional$AQI[df_regional$treated==0 & df_regional$ccc==1])
mean(df_regional$BDR[df_regional$treated==0 & df_regional$ccc==1])
mean(df_local$AQI[df_local$treated==0 & df_local$ccc==1])
mean(df_local$BDR[df_local$treated==0 & df_local$ccc==1])

# state
plots_out_state <- did2s_function(df_state, dvs, dv_colors, 
                                ccc_type = "State-led",
                                custom_title_part = "State-led CCCs:")

aqi_att_state<-plots_out_state$ATT$AQI
aqi_att_state

bad_att_state<-plots_out_state$ATT$BDR
bad_att_state

aqi_plot_state<-plots_out_state$Plots$AQI + ggplot2::scale_y_continuous(limits=c(-13,13),labels=scaleFUN)
aqi_plot_state

bad_plot_state<-plots_out_state$Plots$BDR + ggplot2::scale_y_continuous(limits=c(-0.077,0.077),labels=scaleFUN)
bad_plot_state

# regional
plots_out_regional <- did2s_function(df_regional, dvs, dv_colors, 
                              ccc_type = "Regional-led",
                              custom_title_part = "Regional-led CCCs:")

aqi_att_regional<-plots_out_regional$ATT$AQI
aqi_att_regional

bad_att_regional<-plots_out_regional$ATT$BDR
bad_att_regional

aqi_plot_regional<-plots_out_regional$Plots$AQI + ggplot2::scale_y_continuous(limits=c(-13,13),labels=scaleFUN)
aqi_plot_regional

bad_plot_regional<-plots_out_regional$Plots$BDR + ggplot2::scale_y_continuous(limits=c(-0.077,0.077),labels=scaleFUN)
bad_plot_regional


# local
plots_out_local <- did2s_function(df_local, dvs, dv_colors, 
                                ccc_type = "Local-led", 
                                coeff_range = c(-8, 20),
                                custom_title_part = "Local-led CCCs:")

aqi_att_local<-plots_out_local$ATT$AQI
aqi_att_local

bad_att_local<-plots_out_local$ATT$BDR
bad_att_local

aqi_plot_local<-plots_out_local$Plots$AQI + 
  ggplot2::scale_y_continuous(limits=c(-13,13),labels=scaleFUN)
aqi_plot_local

bad_plot_local<-plots_out_local$Plots$BDR + 
  ggplot2::scale_y_continuous(limits=c(-0.077,0.077),labels=scaleFUN)
bad_plot_local

#####Table 4####
etable(aqi_att_state,bad_att_state, 
       aqi_att_regional,bad_att_regional,
       aqi_att_local,bad_att_local,
       fitstat=~n+ar2, digits =3, digits.stats = 3, tex = TRUE)

#####Figure 7####
figure7<-ggarrange(aqi_plot_state, bad_plot_state,
                    aqi_plot_regional, bad_plot_regional,
                    aqi_plot_local, bad_plot_local,
                    ncol = 2, nrow = 3)

figure7<-plot_grid(figure7)
figure7
ggsave(file = "figure7.jpg", width = 12, height = 12, dpi = 666)
ggsave(file = "figure7.eps", width = 12, height = 12, dpi = 666)

####Mechnism Subgroups####
# goal
mean(df_goal$AQI[df_goal$treated==0 & df_goal$ccc==1])
mean(df_goal$BDR[df_goal$treated==0 & df_goal$ccc==1])

plots_out_goal <- did2s_function(df_goal, dvs, dv_colors, 
                               ccc_type = "Goal",
                               custom_title_part = "CCCs with Clear Goals:")

aqi_att_goal<-plots_out_goal$ATT$AQI
aqi_att_goal

bad_att_goal<-plots_out_goal$ATT$BDR
bad_att_goal

#####Table 7: Panel A####
etable(aqi_att_goal,bad_att_goal, fitstat=~n+ar2, tex = TRUE)

# donation
mean(df_resource_donation$AQI[df_resource_donation$treated==0 & df_resource_donation$ccc==1])
mean(df_resource_donation$BDR[df_resource_donation$treated==0 & df_resource_donation$ccc==1])

plots_out_resource_donation <- did2s_function(df_resource_donation, dvs, dv_colors, 
                                            ccc_type = "Resource: Donation",
                                            custom_title_part = "CCCs Supported by Donations:")

aqi_att_resource_donation<-plots_out_resource_donation$ATT$AQI
aqi_att_resource_donation

bad_att_resource_donation<-plots_out_resource_donation$ATT$BDR
bad_att_resource_donation

# funded by DOE
mean(df_resource_fedfund$AQI[df_resource_fedfund$treated==0 & df_resource_fedfund$ccc==1])
mean(df_resource_fedfund$BDR[df_resource_fedfund$treated==0 & df_resource_fedfund$ccc==1])

plots_out_resource_fedfund <- did2s_function(df_resource_fedfund, dvs, dv_colors, 
                                           ccc_type = "Resource: Federal Funded",
                                           coeff_range = c(-13, 20),
                                           custom_title_part = "CCCs Funded by the DOE:")

aqi_att_resource_fedfund<-plots_out_resource_fedfund$ATT$AQI
aqi_att_resource_fedfund

bad_att_resource_fedfund<-plots_out_resource_fedfund$ATT$BDR
bad_att_resource_fedfund

#####Table 7: Panel B####
etable(aqi_att_resource_donation,bad_att_resource_donation,
       aqi_att_resource_fedfund,bad_att_resource_fedfund, 
       fitstat=~n+ar2, tex = TRUE)

# open membership
mean(df_group_open$AQI[df_group_open$treated==0 & df_group_open$ccc==1])
mean(df_group_open$BDR[df_group_open$treated==0 & df_group_open$ccc==1])

plots_out_group_open <- did2s_function(df_group_open, dvs, dv_colors, 
                                     ccc_type = "Open Membership",
                                     custom_title_part = "CCCs with Open Membership:")

aqi_att_group_open<-plots_out_group_open$ATT$AQI
aqi_att_group_open

bad_att_group_open<-plots_out_group_open$ATT$BDR
bad_att_group_open

# clear benefits
mean(df_group_benefit$AQI[df_group_benefit$treated==0 & df_group_benefit$ccc==1])
mean(df_group_benefit$BDR[df_group_benefit$treated==0 & df_group_benefit$ccc==1])

plots_out_group_benefit <- did2s_function(df_group_benefit, dvs, dv_colors, 
                                        ccc_type = "Clarify Benefits",
                                        custom_title_part = "CCCs with Clear Benefits:")

aqi_att_group_benefit<-plots_out_group_benefit$ATT$AQI
aqi_att_group_benefit

bad_att_group_benefit<-plots_out_group_benefit$ATT$BDR
bad_att_group_benefit

# clear responsibility
mean(df_group_resp$AQI[df_group_resp$treated==0 & df_group_resp$ccc==1])
mean(df_group_resp$BDR[df_group_resp$treated==0 & df_group_resp$ccc==1])

plots_out_group_resp <- did2s_function(df_group_resp, dvs, dv_colors, 
                                     ccc_type = "Clarify Responsibility",
                                     custom_title_part = "CCCs with Clear Responsibility:")

aqi_att_group_resp<-plots_out_group_resp$ATT$AQI
aqi_att_group_resp

bad_att_group_resp<-plots_out_group_resp$ATT$BDR
bad_att_group_resp

# Cross-sector decision
mean(df_group_crosssector$AQI[df_group_crosssector$treated==0 & df_group_crosssector$ccc==1])
mean(df_group_crosssector$BDR[df_group_crosssector$treated==0 & df_group_crosssector$ccc==1])

plots_out_group_crosssector <- did2s_function(df_group_crosssector, dvs, dv_colors, 
                                            ccc_type = "Cross-sectoral Decision Making",
                                            custom_title_part = "CCCs with Cross-sector Decision-makers:")

aqi_att_group_crosssector<-plots_out_group_crosssector$ATT$AQI
aqi_att_group_crosssector

bad_att_group_crosssector<-plots_out_group_crosssector$ATT$BDR
bad_att_group_crosssector

#####Table 7: Panel C####
etable(aqi_att_group_open,bad_att_group_open,
       aqi_att_group_benefit,bad_att_group_benefit,
       aqi_att_group_resp,bad_att_group_resp,
       aqi_att_group_crosssector,bad_att_group_crosssector,
       fitstat=~n+ar2, tex = TRUE)

####Appendix A####
# U.S. Air Quality Index from 1990 to 2020
df_map<-df[,c(3,4,13)]
df_map$fips<-sprintf("%05d", df_map$fips)
df_map$fips<-as.character(df_map$fips)

# Subset data for each year
years <- c(1990, 1995, 2000, 2005, 2010, 2015, 2020)
subset_data <- list()

for (y in years) {
  subset_data[[as.character(y)]] <- subset(df_map, year == y)
}

# Define the colors for the gradient
colors <- c("#00e400", "#ff0000", "#7e0023")
# Define the breakpoints for the colors
breaks <- c(0, 100, 500)

# Generate maps using the subset data
generate_map <- function(year, df_map_year, states, colors, breaks, legend = TRUE) {
  counties_air_year <- plot_usmap(data = df_map_year, 
                                  exclude = c("AK", "HI", "DC"),
                                  values = "AQI", 
                                  color = "black",
                                  size = 0.1)
  
  # Create the base plot
  map_air_year_plot <- ggplot() +  
    geom_polygon(data = counties_air_year[[1]], 
                 aes(x = x, y = y, group = group, fill = AQI),
                 color = "black", size = 0.1) +  
    geom_polygon(data = states[[1]], 
                 aes(x = x, y = y, group = group), 
                 color = "white", fill = alpha(0.01)) + 
    labs(x = element_blank(), y = element_blank(), 
         title = paste("Air Quality in", year)) +
    coord_equal() +
    theme_map() +
    theme_bw() +
    scale_fill_gradientn(colors = colors,
                         values = scales::rescale(breaks),
                         breaks = breaks,
                         labels = c("0", "100", "500"),
                         limits = c(0, 500)) +
    theme(plot.title = element_textbox(size = 22, hjust = 0.5, 
                                       margin = margin(t = 5, b = 5)),
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          plot.margin = margin(t = 0, r = 0, b = 0, l = 0),
          legend.margin = margin(t = 0, r = 0, b = 0, l = 0),
          legend.position = "none")
  
  # Conditional to add the legend if requested
  if (legend) {
    map_air_year_plot <- map_air_year_plot +
      guides(fill = guide_colorbar(
        title = "Air Quality Index",
        barwidth = 20,
        barheight = 1,
        title.position = "top",
        title.theme = element_text(size = 22, face = "bold"),
        label.theme = element_text(size = 22),
        direction = "horizontal",
        frame.colour = "black")) +
      theme(legend.position = c(1.5, 0.2),
            legend.direction = "vertical")
  }
  
  return(map_air_year_plot)
}


maps <- list()
for (y in years) {
  # Include the legend only for the year 2020
  include_legend <- y == 2020
  maps[[as.character(y)]] <- generate_map(y, subset_data[[as.character(y)]], states, colors, breaks, legend = include_legend)
}

# Combine the maps
figure.a1 <- grid.arrange(maps[[as.character(years[1])]],
                             maps[[as.character(years[2])]],
                             maps[[as.character(years[3])]],
                             maps[[as.character(years[4])]],
                             maps[[as.character(years[5])]],
                             maps[[as.character(years[6])]],
                             maps[[as.character(years[7])]],
                             ncol = 2, nrow = 4)

figure.a1<-plot_grid(figure.a1)
figure.a1
ggsave(file = "figure_a1.jpg", width = 15, height = 20, dpi = 666)
ggsave(file = "figure_a1.eps", width = 15, height = 20, dpi = 666)

####Appendix C####
# Descriptive Statistics
df_cov <- data.frame(
  orgtype = df$orgtype,
  AQI = df$AQI,
  BDR = df$BDR,
  pop_density = df$pop_density,
  white = df$white,
  freight = df$freight,
  laborforce = df$laborforce,
  annual_temp = df$annual_temp,
  annual_prec = df$annual_prec,
  total_incen_no = df$total_incen_no,
  hydro_generator_per = df$hydro_generator_per,
  hydro_capacity_per = df$hydro_capacity_per,
  coal_generator_per = df$coal_generator_per,
  coal_capacity_per = df$coal_capacity_per,
  Total_Generators = df$Total_Generators,
  Total_Capacity = df$Total_Capacity,
  total_regulation = df$total_regulation,
  violation = df$violation,
  pro_env = df$pro_env
)

df_cov <- df_cov %>%
  dplyr::rename(
    `Air quality index` = AQI,
    `Bad days ratio` = BDR,
    `Population density` = pop_density,
    `White population (ratio)` = white,
    `Motor freight transportation` = freight,
    `Labor force` = laborforce,
    `Annual temperature (°F)` = annual_temp,
    `Annual precipitation (inch)` = annual_prec,
    `Alternative fuels incentives` = total_incen_no,
    `Ratio of hydroelectric power` = hydro_generator_per,
    `Ratio of hydroelectric capacity` = hydro_capacity_per,
    `Ratio of coal power` = coal_generator_per,
    `Ratio of coal capacity` = coal_capacity_per,
    `Total generator` = Total_Generators,
    `Total power capacity` = Total_Capacity,
    `Regulated facility` = total_regulation,
    `Regulation violation` = violation,
    `Pro-environmental nonprofits` = pro_env
  )

df_cov_ccc<-subset(df_cov, orgtype!=10)
df_cov_npo<-subset(df_cov, orgtype==4)
df_cov_state<-subset(df_cov, orgtype==1)
df_cov_local<-subset(df_cov, orgtype==2)
df_cov_regional<-subset(df_cov, orgtype==3)
df_cov_control<-subset(df_cov, orgtype==10)

# Define the function for descriptive statistics
generate_descriptive_stats <- function(df) {
  # Initialize a data frame to store the results
  descriptive_stats <- data.frame(Variable = character(), Observation = integer(), Mean = numeric(),
                                  SD = numeric(), Max = numeric(), Min = numeric(), 
                                  stringsAsFactors = FALSE)
  
  # Iterate over the variables
  for (var in names(df)) {
    if (var != "orgtype") {
      # Extract data for the variable
      var_data <- df[[var]]
      
      # Calculate descriptive statistics
      mean_val <- round(mean(var_data, na.rm = TRUE), 2)
      sd_val <- round(sd(var_data, na.rm = TRUE), 2)
      min_val <- round(min(var_data, na.rm = TRUE), 2)
      max_val <- round(max(var_data, na.rm = TRUE), 2)
      
      # Add to descriptive stats table
      descriptive_stats <- rbind(descriptive_stats, data.frame(Variable = var, 
                                                               Mean = mean_val, SD = sd_val, 
                                                               Min = min_val, Max = max_val))
    }
  }
  
  # Return the descriptive stats table
  return(descriptive_stats)
}

#####Table C.1####
summary_overall <- generate_descriptive_stats(df_cov)
latex_overall <- xtable(summary_overall)
print(latex_overall,include.rownames=FALSE)

#####Table C.2####
summary_ccc <- generate_descriptive_stats(df_cov_ccc)
latex_ccc <- xtable(summary_ccc)
print(latex_ccc,include.rownames=FALSE)

#####Table C.3####
summary_npo <- generate_descriptive_stats(df_cov_npo)
latex_npo <- xtable(summary_npo)
print(latex_npo,include.rownames=FALSE)

#####Table C.4####
summary_state <- generate_descriptive_stats(df_cov_state)
latex_state <- xtable(summary_state)
print(latex_state,include.rownames=FALSE)

#####Table C.5####
summary_regional <- generate_descriptive_stats(df_cov_regional)
latex_regional <- xtable(summary_regional)
print(latex_regional,include.rownames=FALSE)

#####Table C.6####
summary_local <- generate_descriptive_stats(df_cov_local)
latex_local <- xtable(summary_local)
print(latex_local,include.rownames=FALSE)

#####Table C.7####
summary_control <- generate_descriptive_stats(df_cov_control)
latex_control <- xtable(summary_control)
print(latex_control,include.rownames=FALSE)

####Appendix D####
# Placebo Tests
df_placebo<-subset(df, df$treated!=1)
table(df_placebo$rel_time)
df_placebo$rel_time_placebo<-df_placebo$rel_time+3
df_placebo$treated_placebo<-df_placebo$treated
df_placebo$treated_placebo[df_placebo$rel_time_placebo>=0 & df_placebo$rel_time_placebo!=Inf]<-1

df_placebo_npo<-subset(df_placebo,df_placebo$orgtype==4 | df_placebo$orgtype==10)
df_placebo_gov<-subset(df_placebo,df_placebo$orgtype!=4)
df_placebo_state<-subset(df_placebo,df_placebo$orgtype==1 | df_placebo$orgtype==10)
df_placebo_regional<-subset(df_placebo,df_placebo$orgtype==3 | df_placebo$orgtype==10)
df_placebo_local<-subset(df_placebo,df_placebo$orgtype==2 | df_placebo$orgtype==10)

#####Figure D.1####
plots_out_placebo <- did2s_function_placebo(df_placebo, dvs, dv_colors, 
                                          custom_title_part = "All CCCs:")

aqi_plot_placebo<-plots_out_placebo$AQI + ggplot2::scale_y_continuous(limits=c(-12,12),labels=scaleFUN)
aqi_plot_placebo

bad_plot_placebo<-plots_out_placebo$BDR + ggplot2::scale_y_continuous(limits=c(-0.07,0.07),labels=scaleFUN)
bad_plot_placebo

# plot
figure.d1<-ggarrange(aqi_plot_placebo,bad_plot_placebo,
                        ncol = 2, nrow = 1)

figure.d1<-plot_grid(figure.d1)
figure.d1
ggsave(file = "figure_d1.jpg", width = 12, height = 4, dpi = 666)
ggsave(file = "figure_d1.eps", width = 12, height = 4, dpi = 666)

#####Figure D.2####
# npo
plots_out_placebo_npo <- did2s_function_placebo(df_placebo_npo, dvs, dv_colors, 
                                              custom_title_part = "Nonprofit-led CCCs:")

aqi_plot_placebo_npo<-plots_out_placebo_npo$AQI + ggplot2::scale_y_continuous(limits=c(-12,12),labels=scaleFUN)
aqi_plot_placebo_npo

bad_plot_placebo_npo<-plots_out_placebo_npo$BDR + ggplot2::scale_y_continuous(limits=c(-0.07,0.07),labels=scaleFUN)
bad_plot_placebo_npo

# gov
plots_out_placebo_gov <- did2s_function_placebo(df_placebo_gov, dvs, dv_colors, 
                                              custom_title_part = "Government-led CCCs:")

aqi_plot_placebo_gov<-plots_out_placebo_gov$AQI + ggplot2::scale_y_continuous(limits=c(-12,12),labels=scaleFUN)
aqi_plot_placebo_gov

bad_plot_placebo_gov<-plots_out_placebo_gov$BDR + ggplot2::scale_y_continuous(limits=c(-0.07,0.07),labels=scaleFUN)
bad_plot_placebo_gov

# plot
figure.d2<-ggarrange(aqi_plot_placebo_gov,bad_plot_placebo_gov,
                            aqi_plot_placebo_npo,bad_plot_placebo_npo,
                            ncol = 2, nrow = 2)

figure.d2<-plot_grid(figure.d2)
figure.d2
ggsave(file = "figure_d2.jpg", width = 12, height = 8, dpi = 666)
ggsave(file = "figure_d2.eps", width = 12, height = 8, dpi = 666)

#####Figure D.3####
# state
plots_out_placebo_state <- did2s_function_placebo(df_placebo_state, dvs, dv_colors, 
                                                custom_title_part = "State-led CCCs:")

aqi_plot_placebo_state<-plots_out_placebo_state$AQI + ggplot2::scale_y_continuous(limits=c(-12,12),labels=scaleFUN)
aqi_plot_placebo_state

bad_plot_placebo_state<-plots_out_placebo_state$BDR + ggplot2::scale_y_continuous(limits=c(-0.07,0.07),labels=scaleFUN)
bad_plot_placebo_state

# regional
plots_out_placebo_regional <- did2s_function_placebo(df_placebo_regional, dvs, dv_colors, 
                                              custom_title_part = "Regional-led CCCs:")

aqi_plot_placebo_regional<-plots_out_placebo_regional$AQI + ggplot2::scale_y_continuous(limits=c(-12,12),labels=scaleFUN)
aqi_plot_placebo_regional

bad_plot_placebo_regional<-plots_out_placebo_regional$BDR + ggplot2::scale_y_continuous(limits=c(-0.07,0.07),labels=scaleFUN)
bad_plot_placebo_regional

# local
plots_out_placebo_local <- did2s_function_placebo(df_placebo_local, dvs, dv_colors, 
                                                custom_title_part = "Local-led CCCs:",
                                                coeff_range = c(-5, 2))

aqi_plot_placebo_local<-plots_out_placebo_local$AQI + ggplot2::scale_y_continuous(limits=c(-12,12),labels=scaleFUN)
aqi_plot_placebo_local

bad_plot_placebo_local<-plots_out_placebo_local$BDR + ggplot2::scale_y_continuous(limits=c(-0.07,0.07),labels=scaleFUN)
bad_plot_placebo_local

# plot
figure.d3<-ggarrange(aqi_plot_placebo_state,bad_plot_placebo_state,
                            aqi_plot_placebo_regional,bad_plot_placebo_regional,
                            aqi_plot_placebo_local,bad_plot_placebo_local,
                            ncol = 2, nrow = 3)

figure.d3<-plot_grid(figure.d3)
figure.d3
ggsave(file = "figure_d3.jpg", width = 12, height = 12, dpi = 666)
ggsave(file = "figure_d3.eps", width = 12, height = 12, dpi = 666)

####Appendix E####
# Analysis of States with a History of Poor Air Quality
# exclude NM, AZ, OR, WA, CO, UT,NV, ID, WY, MT, ND, SD, ME
# List of FIPS codes to exclude
exclude_states <- c(35, 4, 41, 53, 8, 49, 32, 16, 56, 30, 38, 46, 23)

# Creating a subset excluding these states
df_estern <- subset(df, !state %in% exclude_states)

df_npo_estern<-subset(df_estern,df_estern$orgtype==4 | df_estern$orgtype==10)
df_gov_estern<-subset(df_estern,df_estern$orgtype!=4)
df_state_estern<-subset(df_estern,df_estern$orgtype==1 | df_estern$orgtype==10)
df_regional_estern<-subset(df_estern,df_estern$orgtype==3 | df_estern$orgtype==10)
df_local_estern<-subset(df_estern,df_estern$orgtype==2 | df_estern$orgtype==10)

# baseline mean
mean(df_estern$AQI[df_estern$treated==0 & df_estern$ccc==1])
mean(df_estern$BDR[df_estern$treated==0 & df_estern$ccc==1])

mean(df_npo_estern$AQI[df_npo_estern$treated==0 & df_npo_estern$ccc==1])
mean(df_npo_estern$BDR[df_npo_estern$treated==0 & df_npo_estern$ccc==1])
mean(df_gov_estern$AQI[df_gov_estern$treated==0 & df_gov_estern$ccc==1])
mean(df_gov_estern$BDR[df_gov_estern$treated==0 & df_gov_estern$ccc==1])

mean(df_state_estern$AQI[df_state_estern$treated==0 & df_state_estern$ccc==1])
mean(df_state_estern$BDR[df_state_estern$treated==0 & df_state_estern$ccc==1])
mean(df_regional_estern$AQI[df_regional_estern$treated==0 & df_regional_estern$ccc==1])
mean(df_regional_estern$BDR[df_regional_estern$treated==0 & df_regional_estern$ccc==1])
mean(df_local_estern$AQI[df_local_estern$treated==0 & df_local_estern$ccc==1])
mean(df_local_estern$BDR[df_local_estern$treated==0 & df_local_estern$ccc==1])

# overall
plots_out_estern <- did2s_function(df_estern, dvs, dv_colors, 
                                 custom_title_part = "All CCCs:")

aqi_att_estern<-plots_out_estern$ATT$AQI
aqi_att_estern

bad_att_estern<-plots_out_estern$ATT$BDR
bad_att_estern

#####Table E.1####
etable(aqi_att_estern,bad_att_estern, fitstat=~n+ar2, digits =3, tex = TRUE)

# npo
plots_out_npo_estern <- did2s_function(df_npo_estern, dvs, dv_colors, 
                                     ccc_type = "Nonprofit-led",
                                     custom_title_part = "Nonprofit-led CCCs:")

aqi_att_npo_estern<-plots_out_npo_estern$ATT$AQI
aqi_att_npo_estern

bad_att_npo_estern<-plots_out_npo_estern$ATT$BDR
bad_att_npo_estern

# gov
plots_out_gov_estern <- did2s_function(df_gov_estern, dvs, dv_colors, 
                                     ccc_type = "Government-led",
                                     custom_title_part = "Government-led CCCs:")

aqi_att_gov_estern<-plots_out_gov_estern$ATT$AQI
aqi_att_gov_estern

bad_att_gov_estern<-plots_out_gov_estern$ATT$BDR
bad_att_gov_estern

#####Table E.2####
etable(aqi_att_gov_estern,bad_att_gov_estern, aqi_att_npo_estern,bad_att_npo_estern,
       fitstat=~n+ar2, digits =3, tex = TRUE)

# state
plots_out_state_estern <- did2s_function(df_state_estern, dvs, dv_colors, 
                                       ccc_type = "State-led",
                                       custom_title_part = "State-led CCCs:")

aqi_att_state_estern<-plots_out_state_estern$ATT$AQI
aqi_att_state_estern

bad_att_state_estern<-plots_out_state_estern$ATT$BDR
bad_att_state_estern

# regional
plots_out_regional_estern <- did2s_function(df_regional_estern, dvs, dv_colors, 
                                          ccc_type = "Regional-led",
                                          custom_title_part = "Regional-led CCCs:")

aqi_att_regional_estern<-plots_out_regional_estern$ATT$AQI
aqi_att_regional_estern

bad_att_regional_estern<-plots_out_regional_estern$ATT$BDR
bad_att_regional_estern

# local
plots_out_local_estern <- did2s_function(df_local_estern, dvs, dv_colors, 
                                       ccc_type = "Local-led",
                                       custom_title_part = "Local-led CCCs:")

aqi_att_local_estern<-plots_out_local_estern$ATT$AQI
aqi_att_local_estern

bad_att_local_estern<-plots_out_local_estern$ATT$BDR
bad_att_local_estern

#####Table E.3####
etable(aqi_att_state_estern,bad_att_state_estern, 
       aqi_att_regional_estern,bad_att_regional_estern,
       aqi_att_local_estern,bad_att_local_estern,
       fitstat=~n+ar2, digits =3, digits.stats = 3, tex = TRUE)

####Appendix H####
# Event-Study Plots for Mechanism Subgroups
#####Figure H.1####
aqi_plot_goal<-plots_out_goal$Plots$AQI + ggplot2::scale_y_continuous(limits=c(-10,5),labels=scaleFUN)
aqi_plot_goal

bad_plot_goal<-plots_out_goal$Plots$BDR + ggplot2::scale_y_continuous(limits=c(-0.05,0.025),labels=scaleFUN)
bad_plot_goal

figure.h1<-ggarrange(aqi_plot_goal, bad_plot_goal,
                     ncol = 2, nrow = 1)
figure.h1<-plot_grid(figure.h1)
figure.h1
ggsave(file = "figure_h1.jpg", width = 12, height = 4, dpi = 666)
ggsave(file = "figure_h1.eps", width = 12, height = 4, dpi = 666)

#####Figure H.2####
aqi_plot_resource_donation<-plots_out_resource_donation$Plots$AQI + ggplot2::scale_y_continuous(limits=c(-11,10),labels=scaleFUN)
aqi_plot_resource_donation

bad_plot_resource_donation<-plots_out_resource_donation$Plots$BDR + ggplot2::scale_y_continuous(limits=c(-0.08,0.07),labels=scaleFUN)
bad_plot_resource_donation

aqi_plot_resource_fedfund<-plots_out_resource_fedfund$Plots$AQI + ggplot2::scale_y_continuous(limits=c(-11,10),labels=scaleFUN)
aqi_plot_resource_fedfund

bad_plot_resource_fedfund<-plots_out_resource_fedfund$Plots$BDR + ggplot2::scale_y_continuous(limits=c(-0.08,0.07),labels=scaleFUN)
bad_plot_resource_fedfund

figure.h2<-ggarrange(aqi_plot_resource_donation, bad_plot_resource_donation,
                         aqi_plot_resource_fedfund, bad_plot_resource_fedfund,
                         ncol = 2, nrow = 2)
figure.h2<-plot_grid(figure.h2)
figure.h2
ggsave(file = "figure_h2.jpg", width = 12, height = 8, dpi = 666)
ggsave(file = "figure_h2.eps", width = 12, height = 8, dpi = 666)

#####Figure H.3####
aqi_plot_group_open<-plots_out_group_open$Plots$AQI + ggplot2::scale_y_continuous(limits=c(-8,8),labels=scaleFUN)
aqi_plot_group_open

bad_plot_group_open<-plots_out_group_open$Plots$BDR + ggplot2::scale_y_continuous(limits=c(-0.05,0.05),labels=scaleFUN)
bad_plot_group_open

aqi_plot_group_benefit<-plots_out_group_benefit$Plots$AQI + ggplot2::scale_y_continuous(limits=c(-8,8),labels=scaleFUN)
aqi_plot_group_benefit

bad_plot_group_benefit<-plots_out_group_benefit$Plots$BDR + ggplot2::scale_y_continuous(limits=c(-0.05,0.05),labels=scaleFUN)
bad_plot_group_benefit

aqi_plot_group_resp<-plots_out_group_resp$Plots$AQI + ggplot2::scale_y_continuous(limits=c(-8,8),labels=scaleFUN)
aqi_plot_group_resp

bad_plot_group_resp<-plots_out_group_resp$Plots$BDR + ggplot2::scale_y_continuous(limits=c(-0.05,0.05),labels=scaleFUN)
bad_plot_group_resp

aqi_plot_group_crosssector<-plots_out_group_crosssector$Plots$AQI + ggplot2::scale_y_continuous(limits=c(-8,8),labels=scaleFUN)
aqi_plot_group_crosssector

bad_plot_group_crosssector<-plots_out_group_crosssector$Plots$BDR + ggplot2::scale_y_continuous(limits=c(-0.05,0.05),labels=scaleFUN)
bad_plot_group_crosssector

figure.h3<-ggarrange(aqi_plot_group_open, bad_plot_group_open,
                      aqi_plot_group_benefit, bad_plot_group_benefit,
                      aqi_plot_group_resp, bad_plot_group_resp,
                      aqi_plot_group_crosssector, bad_plot_group_crosssector,
                      ncol = 2, nrow = 4)
figure.h3<-plot_grid(figure.h3)
figure.h3
ggsave(file = "figure_h3.jpg", width = 12, height = 16, dpi = 666)
ggsave(file = "figure_h3.eps", width = 12, height = 16, dpi = 666)
